//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   CartMesh.cc
 * \author Jeremy Roberts
 * \date   Jul 7, 2011
 * \brief  Member definitions for \ref Mesh
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev:: 141                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-14 19:28:23 +0000 (Wed, 14 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <vector>
#include <cmath>
#include <cstdio>

#include "tran/CartMesh.hh"

namespace slabtran
{

using std::cout;
using std::endl;

// Empty definition for pure virtual base class destructor
Base_Mesh::~Base_Mesh() {}

//---------------------------------------------------------------------------//
// CONSTRUCTORS
//---------------------------------------------------------------------------//

// Constructor using coarse edges and fine counts
template <class D>
Mesh<D>::Mesh(Cell_Edges   &coarse_edges,
              Count_Vector &fine_per_coarse)
  :   Base_Mesh( CountAllCells(fine_per_coarse) )
    , d_dimension(D::dimension)
    , d_number_cells_direction( CountDimCells(fine_per_coarse) )
    , d_region_index(d_number_cells)
{
    int check_dim = 1;
    for ( int dim = 0; dim < d_dimension; dim++ )
        check_dim *= d_number_cells_direction[dim];
    Require ( check_dim == d_number_cells );

    // resize edges, width, and inverse width
    for ( int dim = 0; dim < d_dimension; dim++ )
    {
        Require( coarse_edges[dim].size() ==
                 fine_per_coarse[dim].size()+1 );
        d_edges[dim].resize(d_number_cells_direction[dim]+1);
        d_width[dim].resize(d_number_cells_direction[dim]);
        d_inv_width[dim].resize(d_number_cells_direction[dim]);
    }

    // fine mesh edges and widths
    for (int dim = 0; dim < d_dimension; dim++)
    {
        int fi_store = 0;
        d_edges[dim][0] = coarse_edges[dim][0];
        for (int ci = 0; ci < coarse_edges[dim].size()-1; ci++)
        {
            double w = (coarse_edges[dim][ci+1]-coarse_edges[dim][ci])
                       / fine_per_coarse[dim][ci];
            for ( int fi = fi_store;
                  fi < fine_per_coarse[dim][ci]+fi_store; fi++ )
            {
                d_width[dim][fi] = w;
                d_inv_width[dim][fi] = 1.0/w;
                d_edges[dim][fi+1] = d_edges[dim][fi] + w;
            }
            fi_store = fi_store + fine_per_coarse[dim][ci];
        }
    }

    // generate the region index
    Mesh<D>::GenerateRegionIndex(fine_per_coarse);

}

// Constructor using coarse edges and fine counts, 1d specialization
template <>
Mesh<_1D>::Mesh( Cell_Edges      &coarse_edges,
                 Count_Vector    &fine_per_coarse    )
  :   Base_Mesh( CountAllCells(fine_per_coarse) )
    , d_dimension(1)
    , d_number_cells_direction( CountDimCells(fine_per_coarse) )
    , d_edges( Mesh<_1D>::Vec_Dbl(d_number_cells+1) )
    , d_width( Mesh<_1D>::Vec_Dbl(d_number_cells) )
    , d_inv_width( Mesh<_1D>::Vec_Dbl(d_number_cells) )
    , d_region_index(d_number_cells)
{
    // resize edges, width, and inverse width
    Require( coarse_edges.size() == fine_per_coarse.size()+1 );
    // fine mesh edges and widths
    int fi_store = 0;
    d_edges[0] = coarse_edges[0];
    for ( int ci = 0; ci < coarse_edges.size()-1; ci++ )
    {
        double w = (coarse_edges[ci+1]-coarse_edges[ci])
                   / fine_per_coarse[ci];
        for ( int fi = fi_store; fi < fine_per_coarse[ci]+fi_store; fi++ )
        {
            d_width[fi] = w;
            d_inv_width[fi] = 1.0/w;
            d_edges[fi+1] = d_edges[fi] + w;
        }
        fi_store = fi_store + fine_per_coarse[ci];
    }

    // generate the region index
    Mesh<_1D>::GenerateRegionIndex(fine_per_coarse);
}

//---------------------------------------------------------------------------//

// Constructor using just edges
template <class D>
Mesh<D>::Mesh(Cell_Edges &edges)
  :   Base_Mesh( CountAllCells(edges) )
    , d_dimension(D::dimension)
    , d_number_cells_direction( CountDimCells(edges) )
    , d_region_index(d_number_cells)
{
    int check_dim = 1;
    for ( int dim = 0; dim < d_dimension; dim++ )
        check_dim *= d_number_cells_direction[dim];

    Require ( check_dim == d_number_cells );

    // fine mesh edges and widths
    for (int dim = 0; dim < d_dimension; dim++)
    {
        d_width[dim].resize(d_number_cells_direction[dim]);
        d_inv_width[dim].resize(d_number_cells_direction[dim]);
        for (int fi = 0; fi < edges[dim].size()-1; fi++)
        {
        	double w = (edges[dim][fi+1]-edges[dim][fi]);
        	d_width[dim][fi] = w;
        	d_inv_width[dim][fi] = 1.0/w;
        }
    }
    // generate the region index
    for (int i = 0; i < d_number_cells; i++)
    	d_region_index[i]=i;

}

// Constructor using just edges
template <>
Mesh<_1D>::Mesh(Cell_Edges &edges)
  :   Base_Mesh( CountAllCells(edges) )
    , d_dimension(1)
    , d_number_cells_direction(d_number_cells)
    , d_region_index(d_number_cells)
{
    // fine mesh edges and widths
    for ( int fi = 0; fi < edges.size()-1; fi++ )
    {
    	double w = (edges[fi+1]-edges[fi]);
    	d_width[fi] = w;
    	d_inv_width[fi] = 1.0/w;
    	d_edges[fi+1] = d_edges[fi] + w;
    	d_region_index[fi]=fi;
    }
}

//---------------------------------------------------------------------------//

// Constructor using one width and count
template <class D>
Mesh<D>::Mesh( double    mesh_width,
               int       fine_count )
  :   Base_Mesh( std::pow(fine_count,D::dimension) )
    , d_dimension(D::dimension)
    , d_number_cells_direction(fine_count)
    , d_region_index(d_number_cells,0) // just one region
{
    // resize edges, width, and inverse width
    for ( int dim = 0; dim < d_dimension; dim++ )
    {
        d_edges[dim].resize(d_number_cells_direction[dim]+1);
        d_width[dim].resize(d_number_cells_direction[dim]);
        d_inv_width[dim].resize(d_number_cells_direction[dim]);
    }
    // fine mesh edges and widths
    double w = mesh_width/fine_count;
    double iw = 1.0/w;
    for ( int dim = 0; dim < d_dimension; dim++ )
    {
        d_edges[dim][0] = 0.0;
        d_width[dim].assign(d_number_cells_direction[dim],w);
        d_inv_width[dim].assign(d_number_cells_direction[dim],iw);
        for ( int fi = 0; fi < d_number_cells_direction[dim]; fi++ )
            d_edges[dim][fi+1] = d_edges[dim][fi] + w;
    }
}

// Constructor using one width and count, 1d specialization.
template <>
Mesh<_1D>::Mesh( double mesh_width,
                 int    fine_count )
  :   Base_Mesh( fine_count )
    , d_dimension(1)
    , d_number_cells_direction(fine_count)
    , d_edges( Mesh<_1D>::Vec_Dbl(d_number_cells+1) )
    , d_width( Mesh<_1D>::Vec_Dbl(d_number_cells) )
    , d_inv_width( Mesh<_1D>::Vec_Dbl(d_number_cells) )
    , d_region_index(d_number_cells,0) // just one region
{
    double w = mesh_width/fine_count;
    double iw = 1.0/w;
    d_edges[0] = 0.0;
    d_width.assign(d_number_cells,w);
    d_inv_width.assign(d_number_cells,iw);
    for ( int fi = 0; fi < d_number_cells; fi++ )
        d_edges[fi+1] = d_edges[fi] + w;
}

//---------------------------------------------------------------------------//
// DEFINITIONS
//---------------------------------------------------------------------------//

template <>
const Mesh<_1D>::Vec_Dbl& Mesh<_1D>::edges(int d) const
{
  return d_edges;
}

// Nicely? print me.  Just one long string per dimension.
template <class D>
void Mesh<D>::display()
{
    std::printf("printing the mesh...\n");
    for ( int dim = 0; dim < d_dimension; dim++ )
    {
        std::printf(" dimension: %4i \n",dim+1);
        for ( int i = 0; i < d_edges[dim].size(); i++ )
            std::printf (" %8.4f ", d_edges[dim][i] );
        std::printf("\n");
    }
}
template <>
void Mesh<_1D>::display()
{
    std::printf("printing the mesh...\n dimension: 1 \n");
    for ( int i = 0; i < d_edges.size(); i++ )
        std::printf (" %8.4f ", d_edges[i] );
    std::printf("\n");
}

// Compute total number of cells from count vector
template <class D>
int Mesh<D>::CountAllCells( typename Mesh<D>::Count_Vector &V )
{
    int product = 1;
    for ( int dim = 0; dim < V.size(); dim++ )
    {
        int sum = 0;
        for ( int i = 0; i < V[dim].size(); i++ )
            sum += V[dim][i];
        product *= sum;
    }
    return product;
}
template <>
int Mesh<_1D>::CountAllCells( Mesh<_1D>::Count_Vector &V )
{
    int sum = 0;
    for ( int i = 0; i < V.size(); i++ )
        sum += V[i];
    return sum;
}

// Compute total number of cells by counting fine edges
template <class D>
int Mesh<D>::CountAllCells( typename Mesh<D>::Cell_Edges &V )
{
	int product = 1;
    for (int dim = 0; dim < D::dimension; dim++)
    	product *= (V[dim].size()-1);
    return product;
}
template <>
int Mesh<_1D>::CountAllCells( typename Mesh<_1D>::Cell_Edges &V )
{
    return  V.size() - 1;
}

// Compute total number of cells in each dimension from count vector
template <class D>
typename Mesh<D>::Dim_Vector
Mesh<D>::CountDimCells( typename Mesh<D>::Count_Vector &V )
{
	typename Mesh<D>::Dim_Vector col_sums(0);
	for ( int dim = 0; dim < V.size(); dim++ )
		for ( int i = 0; i < V[dim].size(); i++ )
			col_sums[dim] += V[dim][i];
	return col_sums;
}
template <>
typename Mesh<_1D>::Dim_Vector
Mesh<_1D>::CountDimCells( typename Mesh<_1D>::Count_Vector &V )
{
    typename Mesh<_1D>::Dim_Vector col_sums(0);  // Dim_Vector just an int here
    for ( int i = 0; i < V.size(); i++ )
    	col_sums += V[i];
    return col_sums;
}

// Compute total number of cells in each dimension by counting edges
template <class D>
typename Mesh<D>::Dim_Vector
Mesh<D>::CountDimCells( typename Mesh<D>::Cell_Edges &V )
{
    typename Mesh<D>::Dim_Vector col_sums(0);
    for (int dim = 0; dim < D::dimension; dim++)
    	col_sums[dim]=V[dim].size();
    return col_sums;
}
template <>
typename Mesh<_1D>::Dim_Vector
Mesh<_1D>::CountDimCells( typename Mesh<_1D>::Cell_Edges &V )
{
    return V.size()-1;
}

// Assign region, 3d (default)
template <class D>
void Mesh<D>::GenerateRegionIndex( typename Mesh<D>::Count_Vector &fine_per_coarse )
{
	int i_s = 0;
	for ( int ci = 0; ci < fine_per_coarse[I].size(); ci++ )
	{
		int j_s = 0;
		for ( int cj = 0; cj < fine_per_coarse[J].size(); cj++ )
		{
			int k_s = 0;
			for ( int ck = 0; ck < fine_per_coarse[K].size(); ck++ )
			{
				for ( int fi = i_s; fi < fine_per_coarse[I][ci]+i_s; fi++ )
				{
					for ( int fj = j_s; fj < fine_per_coarse[J][cj]+j_s; fj++ )
					{
						for ( int fk = k_s; fk < fine_per_coarse[K][ck]+k_s; fk++ )
						{
							d_region_index[convert(fi,fj,fk)] =
							ci + cj * fine_per_coarse[I].size() +
							ck * fine_per_coarse[I].size() * fine_per_coarse[J].size();
						}
					}
				}
				k_s += fine_per_coarse[K][ck];
			}
			j_s += fine_per_coarse[J][cj];
		}
		i_s += fine_per_coarse[I][ci];
	}
}

// Assign region id, 2d specialization
template <>
void Mesh<_2D>::GenerateRegionIndex( typename Mesh<_2D>::Count_Vector &fine_per_coarse )
{
	int i_s = 0;
	for ( int ci = 0; ci < fine_per_coarse[I].size(); ci++ )
	{
		int j_s = 0;
		for ( int cj = 0; cj < fine_per_coarse[J].size(); cj++ )
		{
			for ( int fi = i_s; fi < fine_per_coarse[I][ci]+i_s; fi++ )
			{
				for ( int fj = j_s; fj < fine_per_coarse[J][cj]+j_s; fj++ )
				{
					d_region_index[convert(fi,fj,0)] =
							ci + cj * fine_per_coarse[I].size();
				}
			}
			j_s += fine_per_coarse[J][cj];
		}
		i_s += fine_per_coarse[I][ci];
	}
}

// Assign region id, 1d specialization
template <>
void Mesh<_1D>::GenerateRegionIndex( typename Mesh<_1D>::Count_Vector &fine_per_coarse )
{
	int i_s = 0;
	for ( int ci = 0; ci < fine_per_coarse.size(); ci++ )
	{
		for ( int fi = i_s; fi < fine_per_coarse[ci]+i_s; fi++ )
			d_region_index[convert(fi,0,0)] = ci;
		i_s += fine_per_coarse[ci];
	}
}


// explicit instantiations
template class Mesh<_1D>;
template class Mesh<_2D>;
template class Mesh<_3D>;


}  // end namespace slabtran


//---------------------------------------------------------------------------//
//              end of CartMesh.cc
//---------------------------------------------------------------------------//
